%% 含有梯度下降的模拟退火法求解方案
% Grad_Simulated_annealing
function result3=Grad_Simulated_annealing(f,lb,ub)
%% 参数初始化
% lb=[-20,-20];
% ub=[20,20];
n=length(lb);% 变量个数
% alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
% f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
time=3;% 进行3次求解
time1=15;% 每次求解进行15次
time2=1e3;% 求解过程迭代最大次数
result3=zeros(time,n+1);% 储存结果
Tm=100;% 退火初始温度即最大温度
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
% x=x0;
x1=x0;
%% 迭代求解
for m=1:time
    x=lb+m/time*(ub-lb);% 迭代初始点
    result=zeros(time1,n+1);% 储存每次求解结果
    tic
    for k=1:time1
        T=100;% 退火温度初始化
        for n1=1:time2
            L1=f(x);
            k_un=-1+2*rand(1,n);% 随机方向
            z_un=k_un/norm(k_un);% 随机方向规范化
            c_u=zeros(1,n);% 该点梯度初始化
            for bianliang=1:n% 求梯度
                c_un1=x(bianliang)-dx(bianliang);
                c_un2=x(bianliang)+dx(bianliang);
                if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
                    % 量下限时处理
                    r_un=rand();
                    c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
                elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
                    % 于变量下限时处理
                    r_un=rand();
                    c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
                end
                xp=x;x2=x;
                xp(bianliang)=c_un1;x2(bianliang)=c_un2;
                c_u(bianliang)=(f(x2)-f(xp))/(2*dx(bianliang));
            end
            c_u=c_u/norm(c_u+eps);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
            % 梯度方向，否则为次梯度方向。
            randchoice=rand();
            if(randchoice<0.2)
                z_un=-c_u;% 此处说明方向有40%概率为负梯度方向，否则为随机方向
            end
            for i=1:n
                k1_un=x(i)+T/Tm*(ub(i)-lb(i))*z_un(i);
                if k1_un>=lb(i)&&k1_un<=ub(i)
                    x1(i)=k1_un;
                elseif k1_un<lb(i)
                    r_un=rand();
                    x1(i)=r_un*lb(i)+(1-r_un)*x(i);
                else
                    r_un=rand();
                    x1(i)=r_un*ub(i)+(1-r_un)*x(i);
                end
            end
            L2=f(x1);
            if L2<L1
                x=x1;
            else
                r_un=rand();
                p_un=exp(-abs(L2-L1)/T);
                if r_un<p_un
                    x=x1;
                end
            end
            T=(0.99*T);
        end
        if L2<L1
            result(k,:)=[L2,x1];
        else
            result(k,:)=[L1,x];
        end
    end
    result1=min(result(:,1));
    for m1=1:time1
        if result(m1,1)==result1
            result3(m,:)=result(m1,:);
        end
    end
    toc
end